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O As they keep cooling and contracting, Solar System giant planets radiate more en- 

^ ergy than they receive from the Sun. Applying the first and second principles of 

ci thermodynamics, one can determine their cooling rate, luminosity, and temperature 

at a given age. Measurements of Saturn's infrared intrinsic luminosity, however, re- 
veal that this planet is significantly brighter than predicted for its age^ ^. This excess 
luminosity is usually attributed to the immiscibility of helium in the hydrogen-rich 



> 

00 envelope, leading to "rains" of helium- rich droplets^ '*. Existing evolution calcula- 

tions, however, suggest that the energy released by this sedimentation process may 
not be sufficient to resolve the puzzle^. Here, we demonstrate using planetary evolu- 



(^ tion models that the presence of layered convection in Saturn's interior, generated, 

CO like in some parts of Earth oceans, by the presence of a compositional gradient, sig- 

T1 nificantly reduces its cooling. It can explain the planet's present luminosity for a 

. ^ wide range of configurations without invoking any additional source of energy. This 

^ suggests a revision of the conventional homogeneous adiabatic interior paradigm for 

^ giant planets, and questions our ability to assess their heavy element content. This 

reinforces the possibility for layered convection to help explaining the anomalously 

large observed radii of extrasolar giant planets. 

Many arguments suggest the existence of compositional gradients in giant planet interiors, as a 
consequence of either their formation process or their cooling history^^'^^'^^'^-^. However, the effect of 
this gradient on the thermal evolution of the planet is usually neglected for sake of simplicity. When a 
vertical gradient of heavy elements is present in a fluid, the resulting mean molecular weight gradient 
(decreasing upward) can prevent large scale convection to develop by counteracting the destabiliz- 
ing effect of the temperature gradient. The complex interaction between advection and diffusion of 

1 



heat and atomic concentration can trigger a hydrodynamical instability, called the double diffusive 
instability^^, leading to a regime of double diffusive convection (also called semi convection). This 
process significantly affects heat and element transport, as observed in Earth oceans. The instability 
leads to either a state of homogeneous double diffusive convection where diffusive transport is only 
modestly enhanced by small scale turbulence, or to a state of "layered convection", with numerous, 
small convective layers separated by thin, diffusive interfaces corresponding to discontinuities in the 
composition of the fluid, in which transport is enhanced more significantly^^'^^. 

Numerical calculations show that, for relevant thermal and atomic diffusivities, both scenarios are 
possible, depending on the ratio of the compositional to superadiabatic temperature gradients^^. In 
planets, however, this latter quantity is not a free parameter but is imposed by the (measured) energy 
flux to be transported in the planet. It can be shown that this flux is too high to be transported by 
the diffusive regime of the double diffusive instability (said differently, the thermal gradient needed 
to transport the internal flux over the whole planet would be too high to be stabilized by the heavy 
element gradient^-^). Under planetary conditions, the system will thus most likely settle in the regime 
of layered convection. 

To study the impact of layered convection on giant planet thermal structure, we have developed 
an analytical model of the heat transport in such a medium^^. Extending the usual mixing length 
formalism^^ to layered convection, we can calculate the internal temperature gradient as a function 
of the local properties of the fluid and of a unique free parameter, namely the characteristic size of 
the convective layers, /, or equivalently the corresponding dimensionless mixing length parameter, 
a = l/Hp, where Hp is the pressure scale height in the fluid (see Methods). 
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Figure 1: Evolution of the total internal energy (£'tot = fig + ^th + £^rot) with time for various Saturn 
models. The solid curve corresponds to the reference adiabatic model. The long dashed, dashed and 
dotted curves are models with layered convection with a = 10^^'^, 10^-^ and lO^-^-^, respectively. 
Because of the reduced intrinsic luminosity caused by layered convection, these models cool more 
slowly and keep the memory of their initial energetic state much longer. 



A first confirmation of tlie viability of this scenario for giant planet interiors has been given by 
the fact that structure models with layered convection matching all the observational mechanical 
constraints of Jupiter and Saturn (radius, gravitational moments), as well as the atmospheric helium 
and mean heavy element abundances, have been obtained for a rather wide range of mixing length 
parameters, namely a G [10^^ — 4 x 10~^] for Saturn and a G [10^^ — 3 x 10^^] for Jupiter^-^. These 
constraints suggest that layered convection is favored over homogeneous double diffusive convection 
(Supplementary Information). However, the question remains whether layered convection can also 
explain the present luminosity of the giant planets. 

In order to answer this question, we have computed the thermal evolution of planets with layered 
convection. Because departure from adiabaticity can be significant in these interiors, usual isen- 
tropic evolution calculations^ cannot be used. Instead, the evolution is computed by integrating the 
intrinsic luminosity, Ljnt = — d^tot/d^ where the total energy ^tot = ^g H-^th +^rot includes the grav- 
itational, thermal and rotational energies (see Methods). The size of the convective/diffusive layers 
is assumed to have reached an equilibrium value^-^'^^ and is thus kept constant. For all the models 
(i.e. the reference, homogeneous/adiabatic case and the semi convective ones, for any given a), the 
amount of heavy elements is kept constant and equal to the one that fulfills the gravitational moment 
constraints^^ (Supplementary Information). 
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Figure 2: Cooling sequences of Saturn models with layered convection, a) Effective (T^ff; solid 
curves) and intrinsic (Tint; dashed curves) temperature (see Methods) evolution in time for the adia- 
batic reference model (black) and three models with layered convection (a = 10^^'^, 10^^ and 10^^-^ 
from dark to light red), b) Zoom on present era. Dots show the observed effective temperature. At 
early ages the effective temperature of models with layered convection is lower due to inefficient con- 
vection. After a few hundred million years, these models become brighter due to the release of the 
excess of energy stored from the initial state. 



As illustrated in Fig. 1, our results show that starting from the same total internal energy, the time 



required to release a given amount of this energy is significantly longer in the layered case than in the 
homogeneous adiabatic one. This does not imply, however, that the intrinsic luminosity of the object 
will necessarily be larger at any given time. On the contrary, at early ages, objects with a layered 
convection zone are far less luminous because of the reduced heat flux release. However, after some 
time, the decrease in luminosity imposed by layered convection is overpowered by the increase of 
internal energy to be released and planets with layered convection eventually become more luminous 
than the ones with adiabatic interiors (Fig. 2). In Saturn's case, this "crossing time" occurs before 
Saturn's present age so that layered convection yields a larger luminosity at 4.6 Gyr. As shown 
in Fig. 2 and Fig. 3, if the size of the convective layers is about 2-3x10^^ pressure scale heights, 
this process properly leads to the currently observed effective temperature, radius, and gravitational 
moments of the planet without any additional energy source such as helium rains. The radius evolution 
of these models is shown in Supplementary Fig. 5. For this range of layer sizes, the present interior 
temperature of Saturn is only modestly affected compared to the adiabatic case (see Supplementary 
Fig. 6). The temperature during the early evolution, however, is much higher. 
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Figure 3: Impact of the size of the layered convection zone on Saturn's cooling sequence, a) Effective 
(Tgff; solid curves) and intrinsic {T[^u dashed curves) temperature evolution for the adiabatic model 
(black curve) and three scenarios with layered convection, b) Zoom on the present era. From darker 
to lighter red curves we have, i) our baseline scenario with layered convection present throughout the 
gaseous envelope of the planet (Mlc = 0.8A/; a = 2 x 10^^), ii) a scenario where only the inner 50% 
of the envelope exhibits layered convection (Mlc ~ 0.5M; a = 2 x 10^^), iii) a case with a layered 
convection sheU between 0.5 and 0.7M (Mlc ~ 0.2M; a = 10"^). 



Fig. 2 suggests that Saturn models with convective/diffusive layers smaller than 2-3 x 10^ pres- 
sure scale heights will be too bright at the age of the Solar System. This is not necessarily true. 
If layered convection develops only within a restricted fraction of the planet, models with a lower 
convective efficiency (lower a) in the layered zone and with an efficient convection everywhere else 



can also reproduce Saturn's proper cooling timescale. Our results are valid for very different sizes 
of the layered convection zone. Without covering the whole parameter space, we illustrate this by 
showing an evolution track for Saturn where layered convection is restricted to the inner 50% in mass 
of the planet above the core, with a = 2 x 10^^, and an other track where layered convection oc- 
curs in a shell between 50% and 70% of the planet's mass (Fig. 3). This latter case shows that the 
layered convection region does not need to extend down to the core and could be present around the 
molecular/metallic transition region or near an immiscibility region in the gaseous envelope. 

We have also tested the sensitivity of the model to the initial conditions. Indeed, scenarios of 
giant planet formation by core accretion do not provide so far a clear prescription for these initial 
conditions. Large uncertainties in various parameters (gas-to-dust ratio, opacities,...) and the lack of a 
proper treatment of the radiative shock during the gas runaway accretion phase, for instance, hamper 
a proper determination of the planet's initial energy content ^^'^''. Whereas these (unknown) initial 
conditions become inconsequential after about a few tens of million years for an adiabatic planet of 
Saturn's mass^*^, it takes significantly longer (up to a few billion years) for a planet where convection 
is inefficient, because of the lower luminosity and thus the longer thermal, so-called Kelvin- Helmholtz 
timescale {Tkh °^ 1/^int)- Initial conditions thus have a stronger impact on the cooling timescale of 
planets with layered convection. 

As shown in Fig. 4, however, these different initial conditions do not qualitatively affect our re- 
sults, as there is always a possible trade-off between the initial energy content of the nascent planet (at 
the end of the runaway accretion shock) and the thickness or number of convective-diffusive layers. 
As seen in the figure, evolution of models with layered convection can adequately reproduce Saturn's 
presently observed luminosity with very low initial energy content provided sufficiently small convec- 
tive/diffusive layers. Given the current uncertainties of formation models, thermal evolution of Saturn 
with layered, inefficient convection can thus meet the observational constraints for a wide variety of 
initial conditions. Interestingly enough, the same layered convection mechanism, which inhibits the 
heat flow, also provides a plausible explanation to the observed low luminosity of Uranus ^^ In this 
case, however, the "crossing time" discussed above should be larger than the age of the Solar System. 

Because the present luminosity of Jupiter is smaller than the one that is predicted by homogeneous, 
adiabatic evolution models^ (see Supplementary Fig. 7), observations do not seem to support the 
presence of layered convection over the whole planet at the present epoch. This does not preclude, 
however, the possibility for this process either to be present in part of the planet or to have occurred 
in the past, but the planet's gaseous envelope to have been homogenized early enough, due to the 
more vigorous convective flux. The flux emitted by Jupiter during its cooling history is large enough 
to redistribute a significant fraction of a massive initial core (of up to 15 — 20M®) within less than 
4.6 Gyr^^ Layered convection in this planet could then have stopped after the erosion of the whole 
core, or at least of its soluble components'^ '^■'. This would be consistent with both i) the large initial 
core mass needed to form the giant planet on the right timescale"^ and ii) Jupiter's present smaller 
^Qj.g2i,25,i3 jupitej- and Saturn may thus have experienced quite different cooling histories and end 
up with dissimilar interiors. 

Layered convection does not preclude either the possibility of a phase separation between hydro- 
gen and either volatiles or helium in Saturn's interior. In fact, a phase separation favors the occurrence 
of layered convection'^. Indeed, the strong compositional gradients (discontinuities) due to the pres- 
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Figure 4: Past luminosity of Saturn's model with layered convection. Evolutionary tracks of the 
intrinsic luminosity (in units of the presently measured luminosity) of Saturn models integrated back- 
ward in time from presently observed conditions with different mixing length parameters (cu). All the 
models with a > 2 x 10^-^ have too short cooling timescales and cannot be integrated backward for 
4.6 Gyr. Even for quite low initial energy contents, i.e. a low initial intrinsic luminosity, models with 
a sufficiently inefficient convection can reproduce the currently observed luminosity of Saturn. 



ence of an immiscibility region provide very favorable conditions for the development of the double 
diffusive instability and thus of layered convection. Our present treatment of layered convection does 
not take into account the energetic contribution of a possible phase separation. In that case, however, 
the energy release due to the phase separation does not need to be important, as the main role of 
immiscibility would be to trigger the double diffusive instability. Indeed, as both layered convection 
and phase separation increase the planet's cooling time, there is always a possible trade off between 
the respective contributions of these two processes to the planet's cooling rate. Note, however, that 
layered convection yields a hotter interior than the adiabatic evolution, favoring miscibility, especially 
in the past (see Supplementary Fig. 6). As a smaller immiscibility region, if any, is predicted for he- 
lium in Jupiter^'^, this scenario would also explain the dichotomy revealed here between Jupiter and 
Saturn. 

The cooling history of our Solar System giants, and of extrasolar giant planets in general, might 
thus be more complex than assumed by the conventional, simplistic homogeneous, adiabatic interior 
structure and evolution paradigm. The amazing diversity of the observed properties of extrasolar 
planets - in particular the anomalously large observed radii of hot jupiters^^ that cannot be fully 
explained by the various proposed heating mechanisms^^'^^ - combined with the luminosity anomalies 
of our own Solar System giants, clearly suggest a significant revision of this paradigm and point to a 
broader, more complex picture of solar and extrasolar giant planet structure, composition and thermal 
evolution, with a direct impact on giant planet formation conditions. 



Methods 

Layered convection and thermal gradient 

To study the impact of layered double diffusive convection on the interior thermal structure, we have 
developed an analytical model of the heat transport in such a medium^^, which is briefly outlined 
below. 

The fluid is assumed to be composed of a large number of small convective layers of height Z, 
separated by thinner diffusive interfaces. The transport in the convective layers is described by a 
parametrization similar to the standard mixing length theory of convection, for which the mixing 
length is chosen to be equal to the height of the layers (I). The temperature gradient (Vy = ^^^), 
hence the superadiabaticity in the gas, Vt — '^ad, is given as a function of the flux to be transported, 
of the thermodynamic properties of the medium and of /. 

In the diffusive interfaces, the thermal gradient is equal to the gradient needed to transport the 
whole intrinsic flux by diffusive processes alone (Vj). The thermal size of these diffusive interfaces, 
5t, can be estimated by the fact that the convective timescale, given by the inverse of the Brunt- 
Vaisala frequency (Tc = A^^^ ) must be equal to the diffusive timescale in the interface (Tj = 5j/kt, 
where Kj is the thermal diffusivity of the medium). Thus, 5t = \Jkt/Nt, and the only remaining 
free parameter is the height of the convective layers, or equivalently the mixing length parameter, 
a = l/Hp, where Hp is the pressure scale height in the fluid. 

Finally, the mean temperature gradient (Vy) is a linear combination of the two aforementioned 
temperature gradients weighted by the relative size of the convective and diffusive regions, and de- 
pends on the size of the convective/diffusive cells^-' (i.e. on a). 

Planetary model sequences and time integration 

Thermal evolution tracks of a planet of a given mass M, composition (symbolically denoted by an 
array, X, encompassing all the information needed to know the abundances of all the materials at all 
depths, including the core), angular momentum (7), and thermal structure (parametrized here by the 
mixing length parameter a) yield relations of the type 

Lint = ^tcR^OsbT^, = ^TtR^OsB {T^ff-T,\) = Lint(M,l, a, 7,0, 

Em = -Eth +Eg + E^ot = Eioi{M,X, a, 7, t) , 

R = R{M,X,aJj), (1) 

where T^^^ and Tgff are the intrinsic and total effective temperatures, Osb the Stefan-Boltzmann con- 
stant and R the radius of the planet. The mixing length parameter is assumed constant as convec- 
tive/diffusive layers are assumed to have reached a state of equilibrium^^ '^^. Without any further 
guidance from tridimensional hydrodynamical simulations covering a large enough spatial and tem- 
poral domain, it is sensible to start with the most simple assumption. 

The equilibrium temperature (Teq) is the temperature that the planet would have if its intrinsic 



luminosity were zero and represents the contribution to the solar flux through 

a5s7;4=(l-A)L0/(16;ra2), (2) 

where L© is the solar luminosity, a the orbital distance of the planet, and A its Bond albedo. Through- 
out this study, the values used for the equilibrium temperature are T^q = 109 K and 91 K for Jupiter 
and Saturn respectively. This corresponds to an absorbed luminosity of 5.14x lO^^W (1.1 1 x 10^^ W) 
and a Bond albedo of 0.37 (0.32) for Jupiter (Saturn). Within the observational uncertainties, those 
values are consistent with observations. 

The relations described by (1) can be derived with a grid of atmosphere models. These grids of 
atmospheric boundary conditions provide us with the temperature at the 10 bar level (or any reference 
level deeper than the photosphere), rio(7int,^), as a function of the intrinsic effective luminosity (Tint), 
and of the gravity, g = GM/R^. The details of the functional form used to represent the atmospheric 
boundary conditions are detailed in the Supplementary Information. 

Then, for any arbitrary value of the luminosity, Tiq is known, and the T — P profile is integrated 
inward by integrating the standard structure equations for a rotating body (Supplementary Informa- 
tion). Models are computed on a grid of luminosities. For each luminosity, the gravitational, thermal 
{u being the specific internal energy) and rotational energy, 

-^dm, Eth = udm, Erot = i^-r^ (3) 

r J 2 1 

hence the total energy, are computed (see Supplementary Fig. 8). Rotational energy is computed 
assuming negligible angular momentum loss, and is found to represents only a small fraction of the 
total energy (^ 2%). As the ANEOS equation of state provides us with u (the specific internal energy) 
for the volatiles, the thermal contribution of the core is included in our calculations. 

Finally, we are left with a tabulated form of functions of the type Eta = Eio\.{M,X, a,J,L\-at), and 
time dependence can be retrieved by choosing an initial value for the luminosity or initial energy 
content and by integrating 

-Ant = -^- . (4) 
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Supplementary Information 

1 Layered convection or homogeneous double diffusive convec- 
tion? 

The transport efficiency derived from numerical simulations^^' ^^ for a homogeneous double diffusive 
medium (where transport is ensured by small scale turbulence) is equivalent, in our analytical model, 
to a medium in a state of layered convection with a ~ 10^^ — 10^^ (see Figure 3 of Leconte and 
Chabrieri3). 

The fact that measured gravitational moments constrain a to be larger than 10^ {^^) thus seems 
to favor the occurrence of layered convection over homogeneous double diffusive convection in Solar 
System giant planets. In addition, the evolutionary constraint on the luminosity discussed in the 
present study, which constrain a to be larger larger than 10^^ (see Fig. 3), also seems to justify this 
hypothesis. 

2 Atmospheric model grids 

The relations described by (1) can be derived with a grid of atmosphere models. These grids of 
atmospheric boundary conditions provide us with the temperature at the 10 bar level (or any reference 
level deeper than the photosphere), Tio(7int,g), as a function of the intrinsic flux exiting the planet 
(parametrized by the intrinsic effective luminosity, Tint), and of the gravity, g = GM/R^. 

We use the recent atmospheric grids derived specifically for Jupiter and Saturn by Fortney et al.^. 
However, as a planet with layered convection can be quite hot (high internal temperature) but with a 
low luminosity (low effective temperature), our evolutionary tracks can start in the low g, low Tint part 
of the parameter space that is usually not probed by adiabatic evolution models. The grids mentioned 
above thus do not cover this range and extrapolation in g had to be used. The rather low dependence 
of rio(7int,g) on gravity justifies this procedure. 

In order to retrieve the behavior of T\Q{Ti^[,g) in both the high and low Tim regime (arbitrarily 
separated at T^i^, the numerical grids where fitted separately in the two regimes by functions of the 
form 

TxQ{Tm,g)=C + Kg^{T,^,-Toy. (5) 

A smooth function is recovered on the whole temperature range by interpolating linearly between the 
two functions in the range [Tmid — AT, T^i^ + AT] . The 10 bar temperatures derived by this procedure 
do not differ by more than 2-3% with respect to the tabulated values. Values used for the fitting 
parameters are given in Table 1 . 
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Table 1: Parameters used in Eq. (5) to fit the atmospheric grids of ref. 2. 



Jupiter 



Saturn 



Tmid (K) 

AT(K) 



224 
25 



225 
25 





T < Tmid 


T > Tmid 


T < Tmid 


T > Tmid 
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78.4 


-283. 


62.7 


-1680. 
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0.00550 


198. 


0.00202 


831. 


To 


-122. 


143. 


-97.1 


141. 


i8 


-0.182 


-0.114 


-0.180 


-0.0721 


7 


2.09 


0.454 


2.31 


0.293 



3 Planetary structure integration 



When Tio is known, the T 
tions for a rotating body 



P profile is integrated inward by integrating the standard structure equa- 
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dm 
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+ (p(oir), 



Pd 



m 



(6) 
(7) 
(8) 



where m is the mass enclosed in the equipotential of mean radius r, CO the rotation rate of the planet, 
(P(o{r) a second-order correction due to the centrifugal potential. As discussed above, Vj- is either 
equal to Vg^ in the adiabatic case, or to (Vj) which depends on a in a layered convection zone. The 
angular velocity is given by the fact that J = Ico, I being the angular moment of inertia computed for 
each model, is kept constant and equal to the present value. 

A closure equation is provided by the equation of state (EOS) of the mixture along the planet's 
interior profile. Such an EOS is generally given by the so-called ideal volume law for the mixture: 



1 X Y Z 

- = — + — + —, 

P PX PY PZ 



(9) 



where X, Y and Z denote the mass fractions of H, He and heavy elements, respectively. 

For each value of a, the core mass and the mass fraction of heavy elements as a function of 
depth are considered constant in time and equal to the distribution that best matches the gravitational 
moments and observational constraints^^. For adiabatic, homogeneous models, the core masses and 
uniform metal enrichment in the envelope are equal to ~ AM,^ and Z = 0. 11 in Jupiter and 25.5 Mq, 
and Z = 0.05 in Saturn. For models of Saturn with layered convection, the metal enrichments used 
for a given mixing length parameter a are shown in ^^ and the core mass used is ~ 20. SM®. This 
difference in the core mass between the homogeneous model and the layered ones is responsible for 
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the small energy difference remaining at low luminosity in Supplementary Fig. 8. Details about the 
procedure and the equation of states used can be found in ref. 13 and reference therein. 

Supplementary Fig. 7 shows adiabatic, homogeneous evolution tracks for Jupiter and Saturn com- 
puted with this method. As expected, in the homogeneous case, our calculations are in good agree- 
ment with previous results^: models cool to the observed effective temperature in ~ 2.9 Gyr for 
Saturn, and ~ 5.2 Gyr for Jupiter. As we use the same atmospheric model grids, the slight differences 
arising in the cooling times are probably due to the fact that our adiabatic models do not have ex- 
actly the same core mass as the most recent calculations^, and that we take into account the thermal 
contribution of this core. 
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Figure 5: Temporal evolution of the equatorial radius (in units of the measured present equatorial 
radius of Saturn, 60268 km) of an adiabatic Saturn model (solid black curve), and of three models 
with layered convection (a = 3 x 10^^, 10^-^ and 3 x 10^"^ from dark to light red). The effect of the 
deformation due to the fast rotation is taken into account. The time needed for the adiabatic model 
fulfilling the constraints on the measured gravitational moments to cool down and contract to the 
present radius is too short. 
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Figure 6: Internal temperature-pressure profile of our baseline model with layered convection (red 
curves; Mlc ~ 0.8M; a = 2 x 10^-^) at three different ages (dotted: 0.1 Gyr; dashed: 1 Gyr; long 
dashed: 4.6 Gyr). The solid black curve represents the temperature profile inferred for the standard 
adiabatic model. For comparison, the 1 Gyr profile of an adiabatic evolution is shown by the gray 
dashed curve. 



160 

140 
120 h 
100 

80 h 
60 



— *\ ***» — " - 

: ^ i... \ : 

- —--....,- 



2 3 4 

Age (Gyr) 



Figure 7: Effective (T^eff; solid curves) and intrinsic (Tint; dashed curves) temperature as a function of 
the planetary age for our adiabatic reference models of Jupiter (black) and Saturn (gray). Observed 
effective temperatures are also shown as dots. As seen, our calculations in the conventional homo- 
geneous adiabatic case are in good agreement with previous results^: models cool to the observed 
effective temperature in ~ 2.9 Gyr for Saturn, and ~ 5.2 Gyr for Jupiter. 
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Figure 8: Intrinsic luminosity (Lj^t = 4;r7?^Fint) as a function of the total internal energy (Etot = 
Eg + Eth + £'rot) for various sequences of Saturn models. The solid curve is the reference adiabatic 
model. The long dashed, dashed and dotted curves are the a = 10^^'^, 10^^ and 10^^'^ sequences, 
respectively. As expected, for a given internal energy content, the intrinsic luminosity of models 
assuming layered convection decreases when a decreases and is significantly reduced compared with 



the adiabatic case. This increases the Kelvin Helmholtz timescale, Tkh 
needed for the object to loose the memory of its initial energy state. 



GM^/(7?Lint), i.e. the time 
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